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Abstract 



A C++ implementation of the K± jet algorithm for high energy parti- 
cle collisions is presented. The time performance of this implementation is 
comparable to the widely used Fortran implementation. Identical algorith- 
mic functionality is provided, with a clean and intuitive user interface and 
additional recombination schemes. A short description of the algorithm and 
examples of its use are given. 
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1 Introduction 



This paper is intended to be an introduction to the use of KtJet. For a detailed ex- 
planation of the clustering algorithm, its properties, and the physics motivation 
behind it, the reader is referred to [|I| and references therein. The KtJet package 
implements all the features of the Fortran implementation of the i^_L algorithm p[ . 
The philosophy throughout has been to design an interface that users of the Fortran 
code will recognise, but which exploits the advantages of object-oriented design. 
Therefore the names and functions of the input parameters are retained wherever 
possible. The KtJet library, examples of its use and detailed documentation are 
all available on the KtJet website, http://www.ktjet.org/, which also provides 
a link to a CVS repository from which the latest version of the source code may 
be obtained. 



2 The K± algorithm 



The Kj^ algorithm can function in several distinct ways depending on the nature of 
the colliding beams and the physics to be studied. We deal first with hadron-hadron 
and lepton-hadron collisions. There are two modes of operation, the 'inclusive' and 



'exclusive' modes. The difference between the two cases, described in sections 2.1 



and Y,.2[ is in the definition of the hard final state jets, and the separation of these 



jets from the beam remnants. 

For the analysis of e^e~ data, in contrast, there is of course no concept of beam 
remnants. The algorithm proceeds in a very similar way to that employed in a jet 
substructure analysis. We deal with these two cases together in section |2?3|. 



2.1 The inclusive mode 

The algorithm proceeds as follows: 



1. For every final state object[| and for every pair and hi, compute the 
resolution variables dkB and dki- The precise definition of these variables can 
be chosen by the user (using the parameter angle), and will be described in 
detail in section p.4| . They always have the property that in the small angle 
limit they reduce respectively to the squared relative transverse momentum 
of the object with respect to the beam direction, and the squared relative 
transverse momentum of one object with respect to the other. 

^Final state objects could be, for example, partons, hadrons, calorimeter cells, tracks etc. 
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dkB ^ Elelj, (1) 

du ^ mm{ElEf)eli (2) 

At this stage, a dimensionless parameter R is introduced f^, which plays a 
radius-hke role in defining the extent of the jets. This is usually set to 1.0. 

2. Scale the dks by 

dk = dksR^ (3) 

3. Find the smallest value among the dk and dki- If a dki is the smallest, hk 
and hi are combined into a single object with momentum p{^ki) according 
to a user specified recombination scheme (parameter recom) which will be 



described in section p.5[ As an example, recom = 1 would correspond to 
4-vector addition. If a dk is the smallest, object k is defined to be a jet and 
is removed from the list of objects to be merged. 

4. Repeat until all objects have been included in jets. 



2.2 The exclusive mode 

In this mode, the algorithm separates the 'hard final state' from the soft 'beam 
remnants' explicitly. Jets are defined in the hard final state by a stopping parameter 
dcut with dimensions of energy squared. 

1. dkB and dki are defined as in section 

2. Find dmin, the smallest value among the dkB and dki- If dmin > dcut, all 
remaining objects in the event will be classified as jets, and the algorithm is 
complete. 

3. If a dki is the smallest, hk and hi are combined into a single object with mo- 
mentum pi^ki) according to a user specified recombination scheme (parameter 
recom). If a dkB is the smallest, object k is included in a 'beam jet' and 
removed from the list. Go back to (2). 

Note that the stopping parameter d^ut defines the hard scale of the process 
(Aq^^ -C dcut < s, where y/s is the centre-of-mass energy). One can therefore 
think of the hard subprocess being factorised from the low-p± scattering fragments, 
which are thrown away into the beam jets. 
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Instead of specifying the stopping scale dcut, one can choose to stop merging when 
a given number of jets is reached. 



2.3 Subjet analysis and e+e mode 

Once the hadronic final state has been decomposed into jets, the structure of the 
jets themselves can be investigated. This is commonly known as a subjet analysis. 
The procedure is physically (and practically) identical to that employed in the 
analysis of an e~^e~ event. In the e~^e~ case, all final state objects are used as 
input to the algorithm. In a subjet analysis on a particular jet, only the particles 
within that jet are used as input. The steps are as follows; 

1. Define a resolution parameter 

= Ql/El,. (4) 

The parameter E'^^^ may be chosen by the user, but is conventionally taken 
to be the square of the total energy of the e~^e~ event (or the pt of the jet) 
in the frame in which the algorithm is run. These are the default settings in 
KtJet. 

2. For each pair of objects and hi, construct the rescaled resolution variable 
Vki 

yki = dki/El^^ (5) 



where dui is defined as in section 2.1 



3. Find ymim the smallest of the yki- If Vmin < Vcut, and hi are combined into a 
single object with momentum p(fci) according to a user specified recombination 
scheme (parameter re com). 

4. Repeat the procedure until all pairs of objects have yki > ycut- The remaining 
objects are called subjet s. 



As described in section p.2| , one can choose to stop merging when a given number 
of jets is reached, instead of specifying the stopping scale ycut- 



2.4 Jet resolution variables 

The small-angle behaviour of the resolution variables dks and dki given in equations 
m and controls the behaviour of the algorithm in the soft and coUinear limits. 
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Away from the collinear limit, the user has a choice of angular definitions, controlled 
by the parameter angle. Note that the resolution variables are defined such that 
they remain monotonic functions of angle at large angles. Here, and throughout 
this paper, we use rj to denote true rapidity, i.e. rj = 

2.4.1 The angular scheme, angle = 1 



dkB = 2El{l - cosOkB), 

dki = 2mm{ElEf){l-coseki). (6) 

This definition was originally proposed for jets in DIS [^] (where the K± algorithm 
should be implemented in the Breit frame). 



2.4.2 The AR scheme, angle = 2 



dkB - vlki 

dki = mm{plk,pli)Rli, (7) 

where 

Rli = {Vk-ni? + {<Pk-<Pi?. (8) 

This definition corresponds to that used in cone algorithms, and is the most com- 
mon choice for hadron-hadron collisions. 



2.4.3 The QCD emission scheme, angle = 3 

An alternative definition of R\i, motivated by the form of the QCD matrix elements 
for multi-parton emissions: 

Rli = 2[cosh{rik - rji) - cos{(f)k - (pi)]. (9) 



2.5 Recombination schemes 

The recombination scheme defines how two objects hk and hi are merged into a 
single object with 4-momentum p(^ki)- Five schemes are available, selected by the 
parameter re com. 
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2.5.1 The E scheme, recom=l 



Simple 4-vector addition: 

Pki=Pk+Pi- (10) 
This procedure resuhs in massive final state jets. 

2.5.2 The pt scheme, recom=2 



Pt{ki) = Ptk+Pth 

PtkVk + PtiVi 



Vki 



Jkl 



Pt{kl) 

Ptk4>k + pti(t)i 
pt(ki) 



(11^ 



This definition constrains only the 3 spatial components of the object's 4-vector. 
The combined object is made massless by setting its energy equal to the magnitude 
of its 3-momentum. If massive objects are input (for instance particle 4-vectors 
from a simulated event) they are made massless in the same way before the dkB 
and dki are calculated. Compare the Et scheme (section p.5.4| ). 

2.5.3 The scheme, recom=3 

Pt{ki) = Ptk + Pth 

PtkVk + PtiVi 



Tjkl 



■>kl 



2 I 2 ' 
Ptk + Ptl 

Pfk4>k + pl(pi 

pI+pI 



This definition constrains only the 3 spatial components of the object's 4-vector. 
The energy is made equal to to the magnitude of its 3-momentum, thus making 
the combined object massless. Note also that this definition of the scheme is 
that used in the Fortran implementation of the algorithm 0]. It is not equivalent 
to the monotonic pf scheme defined in equations 16 and 17 of reference |jl|. We 
discuss the issue of monotonicity further in section 

2.5.4 The Et scheme, recom=4 

Et(ki) = Etk + Etu 
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EtkVk + Etirji 
Vkl - ^ , 

^t{ki) 

Etk(l)k + Eticpi , , 

(Pkl = r- . (13) 

^t{kl) 

For massless input objects this definition is identical to the pt scheme (section 
|2.5.2| ). It differs solely in the way it deals with massive input objects. If massive 
objects are input, the pt scheme uses their transverse momentum, whereas the Ef 
scheme uses the transverse energy i^sin^. This can have a significant effect for 
low Et jets, particularly for steeply falling distributions close to a cut-off value. All 
combined objects are massless in both cases. 

2.5.5 The scheme, recom=5 



Et{ki) = Etk + Etu 

ElVk + Elr^i 
El + El ' 
%ct>k + Elc^, 



Vkl 



Jkl 



El^k + E^ 
El + E, 



(14) 



This bears the same relationship to the pi scheme that the Et scheme bears to the 
Pt. All combined objects are massless. 



3 Practical implementation of the algorithm 

KtJet uses the HepLorentzVector class of the CLHEP package 0. The 
input 4-vectors (the objects on which the algorithm will be run) may be 
HepLorentzVectors or KtLorentzVectors. The KtLorentzVector class, which 
inherits from HepLorentzVector, carries an internal index which allows the user 
to determine, for example, to which final state jet a particular input particle 
belongs. Examples of the use of KtLorentzVector are given in section |3.2| . 
The output 4-vectors (the jets) are instantiated as KtLorentzVectors. The 
KtLorentzVector class has the constructors 

KtLorentzVector (const HepLorentzVector &) ; 
KtLorentzVector (const KtLorentzVector &) ; 
KtLorentzVector (float Px, float Py, float Pz, float E) ; 
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type 


Beam 


Comments 


1 


ee 




2 


ep 


p in —z direction 


3 


pe 


p in +z direction 


4 


pp 





Table 1: Possible input values for type. 

The first step in running the algorithm is the creation of a KtEvent object. There 
are separate constructors to run the algorithm in inclusive and exclusive modes 
respectively: 

KtEvent (const std: : vector<KtLorentzVector> particles &, int type, 
int angle, int recom, float rparameter) ; 

KtEvent (const std: : vector<KtLorentzVector> particles &, int type, 
int angle , int recom) ; 

There are also versions of these constructors taking vectors of HepLorentzVectors 
instead of KtLorentzVectors. The parameter type should be set to the colliding 
beam type, as defined in table 1. The parameters angle and recom are described 
in sections |2.4| and p.5| . rparameter is the parameter R defined in section p.l| and 
should almost always be set to 1.0, although values smaller than 1.0 have been 
suggested for certain applications 0,0. 

3.1 KtEvent methods 

For the inclusive case, the final state jets are now fully defined, and the jets can 
be recovered using the methods described below. In the exclusive case, the final 
state jets are defined by setting either the stopping parameter dcut- 

void f indJetsD(f loat dcut); 

or by forcing the final state to decompose into jets: 
void f indJetsN(int N) ; 

The dmin value at which the final state changes from + 1 to jets can 
be recovered using the method 
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float getDMerge(int N) const; 

In the case of an e~^e~ analysis, the variable i/cut may be used. The final state jets 
are defined by the method 

void f indJetsY(f loat yCut) ; 

and the i/min value at which the final state changes from + 1 to jets 
can be recovered using the method 

float getYMerge(int N) const; 

By default, i/cut is defined as in equation ^, where Ecut is taken to be the 
total energy in the event. The Ecut value can be set by the user if required: 

void setECut (float ECut) ; 

The KtEvent object has the following methods for recovering its final state jets: 

std: : vector<KtLorentzVector> get Jets () const; 
Returns jets without sorting. 

std: : vector<KtLorentzVector> getJetsEO const; 
Returns jets in order of decreasing energy. 

std: : vector<KtLorentzVector> getJetsEtO const; 

Returns jets in order of decreasing transverse energy. 

std: : vector<KtLorentzVector> getJetsPtO const; 
Returns jets in order of decreasing transverse momentum. 

std: : vector<KtLorentzVector> get JetsRapidityO const; 
Returns jets in order of decreasing rapidity. 

std: : vector<KtLorentzVector> getJetsEtaO const; 
Returns jets in order of decreasing pseudorapidity. 

KtLorentzVector get Jet (const KtLorentzVector &) const; 
Returns the final state jet which contained the input particle KtLorentzVector. 
Note that this method is only available if KtLorentzVectors were used as input 
in the KtEvent constructor. 
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std: :vector<const KtLorentzVector*> & getConstituents () const; 
Returns a vector of pointers to all the input objects in the KtEvent. 

std: : vector<KtLorentzVector> copyConstituents () const; 
Returns a vector of copies of the input objects in the KtEvent. 

int getNConstituentsO const; 

Returns the number of objects in the KtEvent. 

The parameters of a particular KtEvent can be retrieved using the meth- 
ods 

float getETotO const; 
int getTypeO const; 
int getAngleO const; 
int getRecomO const; 
float getECutO const; 
bool isInclusiveO const; 



3.1.1 Monotonicity 

Although in practice most users will not have to face the issue, it is worth bearing 
in mind that recombination schemes 1 to 5 are not guaranteed to lead to mono- 
tonic resolution variables, i.e. it is NOT necessarily true that min{dkB , dki} < 
min{(i'^g, c?^;}, where dk and d'j^ are the resolution variables before and after recom- 
bination respectively. This means that, physically, the question "How many jets 
are there at a particular scale rf?" may not have a unique answer. In KtJet, as 
should be clear from section |2.2| , for a particular value of dcut set by f indJetsD, 
the largest value of N (i.e. the first place at which dmin > dcut) will be returned by 
getNJets. Similarly, there may be no value of dcut for which a particular event has 
jets. There will however always be a dmm value at which the event changed from 
+ 1 to jets, which is the value returned by getDMerge. The above discussion 
also applies to the y variables and associated methods. 

3.2 KtLorentzVector Methods 

These methods would normally be used to investigate the structure and con- 
stituents of final state jets defined in a particular KtEvent. 

const std: :vector<const KtLorentzVector*> & getConstituents () const; 
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Returns a reference to the vector of pointers to all the objects in the 
KtLorentzVector. 

std: : vector<KtLorentzVector> copyConstituents () const; 
Returns a vector of copies of the objects in the KtLorentzVector. 

int getNConstituentsO const; 

Returns the number of objects in the KtLorentzVector. 

bool contains (const KtLorentzVector &) const; 

Check if a jet contains a particular object. For example, 
if (JET. contains (PARTICLE)) { /* do something */ } 

where JET is the KtLorentzVector of a final state jet and PARTICLE is the 
KtLorentzVector of an input object in the KtEvent. 

KtLorentzVector & operator+= (const KtLorentzVector &) ; 

Adds a KtLorentzVector constituent to a jet using the E scheme (4-vector 
addition) and maintains an internal record of constituents (so for example the 
getConstituents method will work on the resulting KtLorentzVector). 

void add (const KtLorentzVector &, int recom) ; 

Adds a KtLorentzVector constituent to a jet using any of the available 
recombination schemes and maintains an internal record of constituents. 

Note that the momentum of a jet (KtLorentzVector) will be given according to 
the recombination scheme used in its construction. If the user wishes to reconstruct 
the momenta according to a different recombination scheme (for example, to re- 
cover the mass for jets which were found using a massless recombination scheme) 
the getConstituents or copyConstituents () methods may be used. The con- 
stituents can then be recombined in a new scheme using the add method. 

4 Subjet analysis 

A subjet analysis is performed on a particular final state jet by constructing a 
new KtEvent object: 

KtEvent (const KtLorentzVector jet &, int type, int angle, int 
recom) ; 

The subsequent analysis and methods arc identical to those described above 
for the e+e" case: for example, subjets can be defined using the findJetsY or 
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f indJetsN methods. E^ut is taken to be the transverse momentum of the jet, 
although this can be changed by the user using the setECut method. 



5 Advanced features of the KtLorentzVector 
class 



The KtLorentzVector class has several features which are used internally by KtJet, 
but which may be useful for the user. 

A KtLorentzVector can be constructed in two different ways. When instan- 
tiated using the constructor KtLorentzVector (const HepLorentzVector 
&) (or KtLorentzVector (float Px, float Py, float Pz, float E)), the 
KtLorentzVector is simply a copy of a single HepLorentzVector with an 
internal index added, i.e. it has no constituents. This index is then used, for 
example, in the getJet method, to ascertain to which final state jet a particular 
input KtLorentzVector belongs. There are two methods which compare the index 
of two KtLorentzVectors: 



bool operator== (const KtLorentzVector &) const; 
bool operator != (const KtLorentzVector &) const; 



These would be useful if, for example, one wanted to ascertain whether the 
highest pt particle in a particular jet was a given input KtLorentzVector. A 
KtLorentzVector may also be instantiated using the constructor 



KtLorentzVector ; 



The KtLorentzVector can then be built up by adding other KtLorentzVectors 
to it using either the assignment operator += or the add method described in 
section |3.2| . The KtLorentzVector will then also carry a list of pointers to the 
constituent KtLorentzVectors. To find out whether a KtLorentzVector has 
constituents, use the method 



bool is Jet const; 



which will be true for a KtLorentzVector with a constituent list. 
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6 Using Kt Jet 



In this section we demonstrate a simple example analysis using the KtJet library. 
We run the inclusive K± algorithm on a proton-proton collision and do a subjet 
analysis on the highest Et jet in the event. 

// This header file must be included. 
#include "KtJet.h" 

std: : vector<KtLorentzVector> jetvec; 

// Loop over input particles, 
for (int i=0; Knpart; i++) { 

KtLorentzVector r = KtLorentzVector (p [i] [0] ,p[i] [1] ,p[i] [2] ,p[i] [3] ) ; 

jetvec .push_back(r) ; 

} 

// We now have a vector of KtLorentzVector, jetvec, 

// containing all the particle 4-vectors, p, in the event. 

// Run the inclusive KT algorithm in PP mode using the covariant E-scheme 
// (type=4, angle=2, recom=l, rparameter=l) . 
KtEvent ev(jetvec,4,2, 1 , 1 . 0) ; 

// Get jets, sorted in Et . 

std: : vector<KtLorentzVector> jetsEt = ev.get JetsEt () ; 

// Perform subject analysis on highest Et jet (angle=2, recom=l) . 
KtEvent jetKjetsEt [0] ,2, 1) ; 

// Decompose jetl into 2 subjets. 
jetl.findJetsN(2) ; 

// Get the KtLorentzVectors of the 2 subjets 

std: : vector<KtLorentzVector> subjetsEt = jetl .get JetsEt () ; 

// Write out the pseudorapidity of the highest Et subjet. 
cout « subjetsEt [0] .etaO « endl; 
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7 Optimisation issues 



The time taken to process events grows with the muhiphcity n of the input objects 
as n^. As the muhiphcity grows with energy, and current and future coUiders 
have ever higher event rates, timing issues are therefore of increasing practical 
importance. 

A few compilation options are worth considering if speed is critical. A compile 
option to switch from single to double precision arithmetic is available, but sin- 
gle precision is used by default as this increases the speed by around 10%. In 
gcc, the option -02 is essential; higher or lower optimisation levels seriously de- 
grade performance. Choice of architecture is also important. For example us- 
ing the gcc compiler flag -march=athlon on an Athlon processor gives a per- 
formance gain of around 5%. More details on performance issues are given at 
http : //www . kt j et . org/. 

8 Adding new functionality 

KtJet interfaces with its recombination and jet resolution schemes via purely ab- 
stract base classes called KtRecom and KtDistance respectively. The packaged 
code comes supplied with 5 Recombination and 3 jet resolution schemes, shown in 
table 2 along with their angle and recom flags. 

Users can define their own recombination and jet resolution scheme classes which 
inherit from the base classes KtRecom and KtDistance and use these instead of 
the supplied schemes. 

To use their own scheme the user needs to pass pointers to the base classes, in- 
stantiated as their own scheme objects, in the KtEvent constructor instead of the 
integer flags. 

As an example, if the user defines a new jet resolution scheme in the class 
KtDistanceNew, and wishes to use the E recombination scheme in an ep collision, 
the inclusive KtEvent constructor should be called as follows: 

KtDistance* distance_scheme = new KtDistanceNew (2) ; 
KtEvent evCparticles , 2, distance_scheme, 1, rparameter) ; 

If the user defines a new recombination scheme in the class KtRecomNew, and wishes 
to use the angle=2 scheme in a pp collision, the inclusive KtEvent constructor 
should be called as follows: 

KtRecom* recom_scheme = new KtRecomNewO ; 
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l\F^^^np 


angle = 1 


KtDistanceAngle(int type) 


angle = 2 


KtDistanceDeltaR(int type) 


angle = 3 


KtDistanceQCD(int type) 


recom = 1 


KtRecomE 


recom = 2 


KtRecomPt 


recom = 3 


KtRecomPt2 


recom = 4 


KtRecomEt 


recom = 5 


KtRecomEt2 



Table 2: The names of the jet resolution and recombination classes included in 
KtJet. 

KtEvent ev(particles , 4, 2, recom_scheme , rparameter) ; 

Users are invited to submit any such additions to the authors for inclusion in future 
releases. 
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